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We study the probability density function P{hm,L) of the maximum relative height hm in a 
^ I wide class of one-dimensional solid-on-solid models of finite size L. For all these lattice models, in 

. the large L limit, a central limit argument shows that, for periodic boundary conditions, P{hm,L) 

takes a universal scaling form P{hm,L) ~ {\/T2w l)~ ^ f{hm/{Vi^iJJL)), with wl the width of the 
fluctuating interface and f{x) the Airy distribution function. For one instance of these models, 
corresponding to the extremely anisotropic Ising model in two dimensions, this result is obtained by 
an exact computation using transfer matrix technique, valid for any L > 0. These arguments and 
exact analytical calculations are supported by numerical simulations, which show in addition that 
the subleading scaling function is also universal, up to a non universal amplitude, and simply given 
by the derivative of the Airy distribution function f (x). 
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I. INTRODUCTION 



The statistics of rare or extreme events in a spatially extended system such as a fluctuating interface has attracted 
considerable recent interest. For example, much theoretical efforts have been devoted to understand the statistics 
r-| , of first-passage events, both in time as well as in space ^ in fluctuating (1 -I- l)-dimensional interfaces of 
Q ' the Edwards- Wilkinson Q or the Kardar-Parisi-Zhang Q varieties. Some of these theoretical predictions have been 
O ' recently verified in experimental systems of monatomic steps on vicinal surfaces 0. Another example involves the 
' study of extremal Fourier intensities for Gaussian interfaces 0- In this paper, we focus on the statistics of yet 

' another extreme observable that has generated considerable recent interest [3, S Il2 Il3 1 namely the statistics of 
J> I the globally maximal relative height (MRH) (relative height refers to the height of the interface measured relative to 
CO ' the spatially averaged height in a finite system) in fluctuating interfaces in its steady state in a finite size system. 

The statistics of MRH in fluctuating interfaces has important theoretical consequences. It turns out that even in 
cases where the steady state of the interface may be simple, such as in a Gaussian interface, analytical calculation of 
the probability distribution of the MRH is nontrivial due to the presence of strong correlations between the heights 
at different points of the interface ^3 ■ While the extreme value statistics of a set of independent or uncorrelated 
\ random variables is relatively easy to study and well understood along with a complete classification of different 
universality classes , much less is known when the random variables are strongly correlated j as in the case of 
■ interfaces. In particular, understanding the issue of universality of the extreme statistics in such a correlated system 
is a challenging problem to which we focus in this paper. 
T ' One of the simplest examples of a fluctuating interface is provided by the (1-1- l)-dimensional Edwards- Wilkinson 
' ^ \ equation Q where the height H{x,t) evolves via a diffusion equation in presence of an additive thermal noise. In 
5 ' the long time limit, the system reaches a steady state which is well known to have the Gibbs-Boltzmann form, 
^ ' Pst oc exp [—Ti] , with a simple Gaussian Hamiltonian 



1 fdH^^ 



where H{x) is the height of the interface and L is the linear size of the substrate. Experimental realizations of 
such a simple Gaussian (1 -I- l)-dimensional interface in thermal equilibrium are provided by thermally fluctuating 
step edges on crystals with attachment/detachment dynamics of adatoms 16]. For periodic boundary condition (pbc) 
H{Q) = H{L), the process H{x) defined by the Hamiltonian in Eq. is just a Brownian bridge. This Hamiltonian has 
a zero mode since it is invariant under H{x) H(x) + c where c is a constant. This is expected since this Hamiltonian 
does not fix the absolute height H of the system. Therefore, a more relevant and physically meaningful observable is 
the relative height, h{x) ~ H{x) — L^^ H[x')dx' , measured with respect to the global spatial average. Although 
the Hamiltonian 7i in Eq. Q retains the same structure in terms of the relative heights, Ti = (1/2) /^^ {dh/dx f' dx, 
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the relative heights /i(a;)'s, by definition, satisfy a global constraint, jj^ h{x)dx =^ 0. This constraint turns out to play 
a crucial role in determining the MRH distribution llfll] . 

The variable h{x)^s at different space points can be shown to be strongly correlated. For example, for the pbc, the 
two point correlation function is given by [10| 



C{x, L) = {h{xQ)h{xo + 2;)) = ^ 
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The onsite variance {h?'{x)) = C{0,L) = L/12 and the width W]^ = y/ {h'^{x)) — y/L/12. Thus the typical relative 
height scales as /i ~ L^/^. Let hm denote the maximum relative height in a given sample, i.e. hm = max[{h{x)} , < 
X < L]. Clearly hm varies from sample to sample. What is its probability density P{hm,L)? Since the relative 
heights h{x) in the steady state are strongly correlated random variables, the calculation of P{hm,L) is nontrivial. 
In Ref. m, it was demonstrated numerically that P{hm,L) has a scaling form, P(/im,L) ~ L~^/'^ f{hm/ VL). Later, 
it was proved that for a continuous interface, this scaling holds for all L (and not just for large L). Expressed in 
terms of the average width wl — yjL/\2, 

for all L. Moreover, the scaling function /(x) was computed exactly 0,^3 using path integral techniques and it was 
shown to be sensitive to the boundary conditions. For example, for the pbc, it was shown that the Laplace transform 
of f{x) is given by ^ 

00 

f{x)e-'''dx = sV2^J2 e-"^-^'^'2-i/^^ (4) 

k=l 

where a^'s are the magnitudes of the zeros of the standard Airy function Ai(z) on the negative real axis [Tt'I. For 
example, ai = 2.3381 . . . , ^2 = 4.0879 . . . , 03 = 5.5205 ... etc Note that the function f{x) defined in Eq. Q 

can be interpreted as a normalized probability density function (pdf) since f(x ) > for all x > and it can be shown 
that f(x)dx = 1. The function f{x) has the following asymptotic tails [i3 

fix) ~ X-5e-2a?/27.^ as X^O 

f{x) e"*^^^ as X ^ 00. (5) 

Remarkably, this function f{x) app eared before in a number of seemingly unrelated problems in computer science 
and graph theory and is known |lCt Il8l Il9j| as the Airy distribution function (not to be confused with the Airy 
function Ai(x) it self) . Essentially the function f{x) describes the pdf of the area under a Brownian excursion on an 
unit interval [23, ^M- Thus it is interesting that the same pdf also describes the MRH distribution for a Gaussian 
interface with pbc^, T^. The results in Eqs. (j^J and Q provide one of the rare exactly solvable cases for the 
distribution of the extremum of a set of strongly correlated random variables. 

A natural question is to what extent this MRH distribution P{hm,L) and the associated scaling function f{x) is 
universal? We have already mentioned above that it is sensitive to the boundary condition. This is natural since hm 
is the global maximum over the full sample. However, suppose we fix the boundary condition to be, say periodic, and 
ask how sensitive P{hm,L) is on the details of the short range interaction in the Hamiltonian Ti, in Eq. in the 
steady state? In this paper we address this issue and study, both analytically and numerically, the MRH distribution 
in a class of one dimensional solid-on-solid (SOS) models The SOS model is defined on a lattice with L sites 

labelled 0, 1, ... L — 1, with pbc (so that the site L is identified with site 0). Note that we have set the lattice constant 
to be unity so that L also represents the total length of the substrate as in the continuous Hamiltonian in Eq. 
The model, at equilibrium, is described by the Hamiltonian 

Hp = KY,\H,-H,+^\P (6) 

where iJ^ e] — 00, -|-oo[ is a continuous height variable and p > is a positive parameter. This model is a continuous- 
height version of the discrete- height SOS models where the height variables HiS are integers For the discrete- 
height SOS model, one can generate a family of models by tuning the parameter p. For example, the case p = 1 
corresponds to the extremely anisotropic Ising model in 2 dimensions, where as the case p — > cx) corresponds to the 
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restricted SOS model where neighbouring heights can differ at most by one unit |22|. Several authors have studied 
the p = 1 model, both for discrete heights ,2^] and for continuous heights in the context of wetting phenomena 
in two dimensions. 

In this paper we restrict ourselves, for simplicity, to the model in Eq. ^ with continuous height i/^'s, though 
our main asymptotic results will be valid even for the discrete-height models. We define the relative heights hi — 
Hi — '^iHi/L and study the distribution of the MRH P{hm,L) in this model. Note that the case p = 2 in Eq. 
© corrresponds to the spatially discretized version of the continuous-space Gaussian interface model in Eq. JQ. 
Therefore, in this case, one would expect that as the number of sites L becomes asymptotically large, one would 
recover the continuum limit results for the MRH in Eqs. lO and Less obvious is what happens when p 0. Our 
main object here is to study the MRH distribution for arbitrary p > 0, to see how it depends on the parameter p and 
also to investigate the leading finite size effects. 

Our main results are summarized as follows. We will demonstrate that asymptotically for large L, the MRH 
distribution P{hm,L) is described by the scaling form in Eq. Q where the width wl depends on p (weakly) but 
the scaling function f{x) is universal, i.e. independent of p and is described by the Airy distribution function in Eq. 

This result is less obvious apriori for p ^ 2. We will first present a general argument based on the central limit 
theorem that will suggest this universality with respect to p. Next we will present numerical simulations in support 
of this universal result. In fact, this central limit argument holds even for a more general class of SOS models defined 
by the Hamiltonian 

exp[-7^] (7) 

i 

where g{x) is an arbitrary positive, symmetric and normalized (to unity) function, but with a finite second moment 

For the special case p = 1 of the model in Eq. ©, we will present an exact calculation of P{h„i, L) for any L that 
proves explicitly this universality in the large L limit. Moreover, for p = 1, we also calculate exactly the subleading 
correction to the leading scaling form for large L and then demonstrate numerically that even the subleading scaling 
function, up to a p dependent amplitude, is also universal with respect to p and the associated subleading scaling 
function is simply the derivative f'{x) of the Airy distribution function in Eq. 

The paper is organized as follows. In the next section (II), we present an argument for the universality of the 
asymptotic MRH distribution based on the central limit theorem. In section III, we present a general set-up to 
compute the MRH distribution in the general SOS models defined in Eq. ^ for arbitrary L. Section IV contains an 
exact result, valid for all L, for the p = 1 case. In section V, we present the details of the numerical simulations. We 
conclude in section VI with a summary and outlook for future studies. Finally, in appendix A, we present an exact 
calculation of the width of the interface defined in Eq. {Tj). 



II. A GENERAL UNIVERSALITY ARGUMENT BASED ON CENTRAL LIMIT THEOREM 



Our starting point is the general Hamiltonian in Eq. ((Tj). The product form on the right hand side of Eq. ((7|) 
indicates that one can interpret the successive height differences -ff^+i — Hi = ^^'s to be a set of random numbers 
drawn from the joint distribution 



(8) 



where is such that the joint distribution is normalized and g{x) is a symmetric, normalized pdf with a finite second 
moment cr^, but otherwise arbitrary. The delta function ensures that the path is periodic, Hq = H^. Therefore, the 
height profile Hi can be interpreted as a discrete random 'bridge' process, bridging the two ends Hq — Hl and in 
between, performing a discrete random walk by choosing at each step a random jump variable ^ drawn independently 



from the normalized pdf Prob(^ 
© the noise pdf is 



x) = g{x) with a finite variance a . For the special case of the Hamiltonian in Eq. 



Prob(C 



gp{x) 



'K\x\' 



where the normalization constant i?„ = f°° dxe ^l^l'' and the variance ct^ — f°° 

P J— oo P J — c 



2T{l/p) 



K-i/p 



r(3/p) 
r(i/p) 



j{x) dx are obtained as 



(9) 



(10) 
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In particular, Bi = 2/K and af = will be used later in section IV. 

In the limit of large L, a standard central limit argument asserts that the discrete random bridge process Hi will 
converge in law, up to a nonuniversal scale factor that depends on the variance of the noise pdf g{x), to the 
continuous Brownian bridge process H(x) described by the Hamiltonian in Eq. Of course, this will be true as 
long as Prob(^ — x) = g{x) decays faster than for large i.e. the second moment of g{x) is finite. Naturally 
the same conclusion will hold for the relative height hi — Hi — '^iHi/L as well. More precisely, for large hi and 
large L, the dimensionless scaled variable hi/w^ (where wl is the width defined via the relation wj^ — '^i{hi) / L) will 
have identical asymptotic statistics as that of its continuous counterpart h(x) j ^ L jVl. Naturally the scaled MRH 
hm/wL where h^ = max ({^i}, i = 0, 1, • • • , L — 1) is expected to follow the same asymptotic statistics as that of its 
continuous counterpart, as long as is finite. This argument thus predicts that for arbitrary symmetric g{x) with a 
finite variance tr^, the asymptotic MRH distribution will converge to 

P{^h^,L)^^^f(^^\ (11) 



l2wL \VT2wlJ 

where f{x) is the Airy distribution function defined in Eq. Q). 

Indeed, 'nonuniversality' of the MRH distribution, i.e. its dependence on the noise pdf Prob(^ = x) = g{x) enters 
only through the width wl in Eq. (|ll|l . In fact, for the discrete process with arbitrary noise pdf g{x), one can 
calculate the width wl exactly (as derived in appendix A) for all L and one gets 



L'-i. . /-°L[g(fc)]^-^ [g'{k)Ydk 

where g{k) — g{x)e^''^ dx is the Fourier transform of g{x) and g'{k) = dg/dk. In particular, for large i, one can 

show (see appendix A) that « cr^/L where = ~ff"(0) is the variance of g{x). Then wl ~ <7y/L/12, i.e. the 
same result as in the continuous case, apart from the nonuniversal scale factor a. 

This completes our argument for the universality of the MRH distribution in the SOS model defined in Eq. {Tj) 
which predicts that, for large L, the MRH distribution for arbitrary symmetric g{x) > (with a finite second moment 
a^) would have the scaling form in Eq. Hll|l with wl ~ a^J L/\2 and with a universal scaling function f{x) (Airy 
distribution function) as defined in Eq. 10}. The universal asymptotic scaling function predicted in this section are 
further supported by the exact result in section IV for the p = 1 case of the model in Eq. ((BJ and numerical simulations 
for other values of p as presented in section V. 



III. THE GENERAL SET-UP TO COMPUTE THE MRH DISTRIBUTION FOR ARBITRARY L 



Let us outline the general set-up to calculate the MRH distribution in the SOS models defined by the Hamiltonian 
in Eq. 10. This set-up is a discrete analogue of the method used for the continuous Gaussian interface model 9, 10]. 
In terms of the relative heights hi = Hi — J^i^i/-^^ ^he SOS Hamiltonian retains the same form as in Eq. ||7J). 
However, the relative height /i^'s now, by definition, satisfy a global constraint, hi = 0. Moreover, the pbc requires 
ho = hL. Thus, the joint probability distribution of the relative heights {hi} with i = 0,1,...,L — 1 in the steady 
state can be written as 



L-1 



lI) 



i=0 



L-1 
.i=0 



(13) 



where the delta function on the right hand side ensures the global constraint hi — 0. The normalization constant 
Al, to be calculated later, is fixed by demanding that the joint distribution P[{hi}] is normalized to unity. 

We wish to compute the pdf P{hm, L) of the MRH h,n — rnax [/iq, /ii, . . . , It turns out to be easier first 

to compute the cumulative distribution F{hm,L) — Prob [max({/ii}) < h„i,L]. The pdf of the MRH is simply the 
derivative, P{hm,L) = dF{hm, L) /dhm- Note that F{hm,L) is the probability that the maximum of all heights are 
less than hm, which is the same as the probability that all the height variables hi < h^. Using the joint probability 
distribution of /i^'s in Eq. H13|l . one can then express F{hrn,, L) as a multiple integral 



F{h^,L) 



dho / dhi , 



dhL-lP[{h,}] 



A, 



Vh 



Y[g{\h^ - h,+i\) 



s(Y^h. 



(14) 
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where Dh is a shorthand notation for Vh = dho dhi . . . d/iL-i- Making a change of variable t/i = hm — hi, Eq. 114() 
becomes 



/•oo 

F{h^,L)^AL / Vy 
Jo 



Y[9{\yi-yi+i\) 



siY^y^-f^rnL] . (15) 



Note that the h„i appears on the right hand side only inside the delta function as a combination A = hj^L. Fur- 
thermore, can not be negative as this would indicate that all the relative heights hi < which can not be true 
due to the exact constraint hi = 0. Thus, the pdf P{hrm L) has nonzero support only for < h^ < oo. Writing 
F{hm,L) = T{A, L) and taking its Laplace transform with respect to A wc get 

/ T{A,L)e-^^dA^ALZ{\L) ; Z{\L)^ Vy {\y, - y,+i\) e-^^^ . (16) 
Jo Jo i 

For the special case when g(x) oc exp[— i4r|a:;|P] as in Eq. 0, the partition function is given by 



Z{\ L) = Vy exp (^K ^ \y, - T " ^ E 



(17) 



Thus if we can determine the partition function Z{\, L), then by inverting the Laplace transform we get F(hm, L). 
The normalization constant A^ can then be fixed by using the fact that F[hm ~^ oo, L) = 1 (which follows from the 
fact that the MRH pdf must be normalized to unity). In fact, this procedure can be directly implemented in Eq. 116|) 
by taking the A ^ limit. As A ^ 0, T{A, L)e-^^ dA 1/A using the fact F{A ^ oo,L) = 1. Then it follows 
that 

Al = lim ^ . (18) 

The partition function Z{X,L) in Eq. (|16|l can, in principle, be computed by a transfer matrix method for a general 
g{x). However, an explicit exact solution for all L can be obtained for the p = 1 case in Eq. as demonstrated in 

the next section. 

We end this section by pointing out an interesting connection between the MRH distribution in Eq. (|14|l and a 
first-passage or no barrier crossing probability of a discrete random acceleration process. It is evident from Eq. H13() 
that the relative height variables hi also form a random bridge (pbc) as the absolute height HiS, except with an 
additional long range constraint = as manifest by the delta function in the joint distribution in Eq. I|13|l. 

Thus the process hi, though locally a random walk, it is conditioned to return to its initial point (pbc) and also 
has to remember that at the end of L steps the total 'area' under the walk must be identically 0. This long range 
global constraint thus induces a 'memory' effect and makes the evolution of hi a non-Markovian process. This fact 
was already noted in the context of the continuous Gaussian interface model 0, 0| . This non-Markovian process can 
be recast as a Markov process via the following trick. We define an additional 'area' variable Ai = X]j=o ^i-i with 
Aq = 0. Then the joint two variable process {hi, Ai) actually evolves by a Markov process 

hi+i = hi +^i , Ai+i = Ai + hi (19) 

where ^^'s are drawn from the joint distribution in Eq. ((SJ. 

Thus, Eq. H19|) is just the discrete version of the random acceleration problem [25l {i.e. the second derivative 
of A is noise). This joint process {hi,Ai) starts with its initial value {h(),Ao = 0) at step i = and ends at its final 
value {h^ = ho,AL = 0) after L steps ensuring both the pbc and the zero area constraint at the end of the process, 
Al — J2i=o ~ ^- This then completely describes the inferface height profile hi at equilibrium. The cumulative 
MRH distribution F{hm,L) in Eq. H14|l is then just the probability that the joint process {hi,Ai) evolving via Eq. 
(|19|l reaches from its initial value (/iq, = 0) to its final value {h^ — ho,AL — 0), but with the restriction that the 
/li's stay below the level hm up to L steps. This is then just a restricted propagator of the discrete random acceleration 
process. 

IV. EXACT SOLUTION FOR ARBITRARY L USING TANSFER MATRIX FOR p = 1. 

The partition function Z{X,L) in Eq. (|16l) can in principle be computed, for arbitrary g{x), using a standard 
transfer matrix technique. The transfer matrix T associated to Z{X, L) in Eq. (|16() is defined via 

{x\T\y) = exp {-Xx/2) g{\x - y\) exp (-Ay/2) (20) 
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and the partition function is simply the trace, Z{X,L) = Trace[r^]. Let us define as its eigenfunctions that 
satisfy the eigenvalue equation 



dy {x\T\y) <j,iy) = Ecj){x) (21) 
where E label the eigenvalues. Making a substitution, (j){x) = tp{x) e'^^l"^ one gets an integral equation 

dyg{\x - y\) exp (-Ay) ^^(y) = Ei,(x) (22) 



For an arbitrary g(x\ the integral equation H22|l is hard to solve. However, one can make progress for the special 
case of the Hamiltonian in Eq. © with p = 1 that corresponds to g{x) = exp[— ii'lxl] to which we now focus. For this 
case, the integral equation H22(l can be transformed into a Schrodinger-like (though not quite the same) differential 
equation by using the identity {—d^/dx'^ + i-f^-jg-xix-yl _ 2Kd{x — y) Q^. We get 



^"i^) + ( ^e-^- - ] ^(x) = (23) 



with the boundary conditions 



(i) lim ^P{x) = , i^^)t^ = K (24) 

which follow from Eq. (|22|l . 

The differential equation H23|) can be further reduced to a standard Bessel form by a change of variable. The general 
solution to this equation has two linearly independent parts, but one of them is not allowed due to the condition (i) in 
Eq. 12311. One finally gets the solution of Eq. ^ as tp{x) = AJ2K/\[{'2/>^)W2K/E)e~'^'-'^'^], where J^, (x) is a Bessel 
function and A is an arbitrary amplitude. The 2nd boundary condition (ii) in Eq. (|24|) determines the eigenvalue E 
as a root of the following equation 

J2K/xiu)/J2K/x{u) = -2K/{\u) (25) 



where u = 2^j2K/E/\. Using the identity Jl{x) — J^-i{x) — {v /x)Ji,{x), Eq. I^K]) becomes: 



2 / A 



J2K/\-l{u) = , i.e. En = -^j2K/\-l,n ) (26) 



where j^^n is the n root of Jv{x) on the real axis. Thus Z{X,L) is given by the exact formula, valid for all L, 

-2L 



Z{\, L) ^B^Y^i WTPJ^K/x-Un ) (27) 



where Bi = 2/K. This completes our derivation of the exact partition fimction Z{X, L) for all L. 

Asymptotic results for large L: We next focus on the large L limit. If one naively takes the large L limit of 
Z{X, L) keeping A fixed, it is just enough to retain only the largest eigenvalue of the transfer matrix in the partition 
sum. However to extract the leading scaling behavior of the cumulative distribution F[hm,L), we need to take the 
limits L CO, hm oo simlutaneously while keeping the ratio hm/VL fixed. Since A = hm L ~ L^^^ is the variable 
conjugate to the Laplace variable A in Z{X, L), the appropriate scaling limit is L — > oo, A — > but keeping the scaling 
combination XL^^^ fixed. Hence we define a new variable 

/2 

s = (TiXL^''^] where = (28) 

K 

We next expand (y^i2K/\~i,n) at large L, keeping s fixed. To proceed, we use the expansion of jv,n for large 
n fixed [13: 

J.,n^^+^>^'^' + Oi^-'/') (29) 
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where a„ is the amphtude of n*"^ zero of the Airy function Ai(z) on the negative real axis. Using this expansion H29|l . 
we obtain 



Z{X, L) = B^Y^ e-'^.s^^^'-'' + ^ + OiL-^)j . (30) 
Substituting this result in Eq. (|17fl . one gets for p — 1 



J T{A, L) e-^^ dA = Al ^ ^-^^s-/-^^'" j^i + + 0(^-1; 



(31) 



so that, using the Bromwich formula for Laplace inversion and substituting A = h„iL one obtains the cumulative 
distribution as 

r-Ao + ioo ^ ^ ^r3/2^2/3^-l/3 '^^^^^^^ + 0{L^^)] (32) 



J\o-roo 27ri ^ V 



2L 



where the integration is along any imaginary axis whose real part Aq must be to the right of all the singularities of 
the integrand. Performing the change of variable s = aiXL^/^ and taking the derivatives with respect to hm one finds 
that P{hm,L) has the structure 



P{h^.L)- ^^^^ 



where f{x) is the leading scaling function in the large L limit. From Eq. (|32|) . its Laplace transform f{s) is given by 

/(s) = / dx fix) e-^^ = ,s^/2^ V e-""^'''2-/^ (34) 

as in Eq. Q and thus f{x) is the Airy-distribution function. Similarly it is easy to check from Eq. (|32|l that 

fM^ dxh{x)e-''- f{s) . (35) 

We are now ready to determine the normalization constant A^. Since the pdf P(h„i,L) is normalized to unity, 
we have P{h,n, L) dh„i — 1. Integrating Eq. and using the facts: (i) f{x)dx ~ /(O) = 1 (which can be 
shown from Eq. and (ii) /i(x) dx — /i(0) = (as follows from Eq. (jSHll), we get 

-AL^[l + OiL-')] = l (36) 



(33) 



that gives, using Bi — 2/K, 



Al = V2^a, L^'^ [-] [l + 0{L-^)\ . (37) 



Substitution in Eq. H33|) then gives the asymptotic scaling function and its first corrrection 

P(,„., . ( J^) + _L/, + .. ,38, 

where f{x) and fi{x) are given by their Laplace transforms in Eqs. (|34|l and (|35|l respectively. This is our main result 
of this section. 

It turns out that the Laplace transform in Eq. H34() can be formally inverted [lol l2l| 

/(^) = ^Ee-''"/^'6^/3t/(-5/6,4/3,5„/x2) (39) 

n—1 

where = 2a^/27 and U{a^b,z) is the confluent hypergeometric function This inversion is useful as it can 

then be evaluated using Mathematica and subsequently one can compare it with the numerical simulations later. 
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The leading 'correction to scaling' function fi{x) has received some attentions recently in a class of random walk 
problems that arise in the context of rectangle packing problem |27|. |28l as well as in the problem of computing the 
total flux of diffusing particles to a spherical trap in 3 dimensions [23|- It is interesting that in the present problem 
we can compute this correction to scaling function fi{x) exactly. Inverting the Laplace transform in Eq. (|35ll . we get 
= f (x) / where f (x) can be evaluated by taking the derivative of Eq. H39|l and using properties of C/(a, 6, z). 
This gives 

/'(^) = -y /(^) + ^ E e-""/^'^'/' (ui-5/6, 4/3, bjx') - ^[7(1/6, 7/3, bjx^)) (40) 

^ n=l ^ ^ 

Finally, using ^/T2wl — aiL^^^ + 0{L^'^^^) from Eq. H12() . we can express the scaling in terms of the width wl 
instead of L^^^. The first two terms are given by 

In particular, from Eq. H41|) . one obtains the mean value E[hm\ as 

E[h^]^a,VL-^+0{L-^/^) (42) 

The leading correction to the asymptotic L^/^ behavior of the expected MRH turns out to be a negative constant, as 
in the case of the expected maximum of an ordinary random walker treated recently in Refs. j27ll28ll2ij'|. 

The exact solution presented in this section proves that indeed the leading scaling function of the MRH pdf in the 
SOS model with p = 1 is given by the Airy distribution function, in agreement with our general argument presented 
in section II. In the next section we present numerical simulations for other values of p to provide further support in 
favour of the universality of the Airy distribution function. In addition, the exact solution for p = 1 also provides the 
leading correction to scaling function /i (x) exactly. We will see in the next section that simulations for other values 
of p indicate that even this leading correction to scaling function /i(x) is also universal, i.e. independent of p. 



V. NUMERICAL SIMULATIONS 



In order to generate numerically height profiles Hi (with pbc) according to the Boltzmann-weight associated to 
Tip in Eq. © we use the algorithm proposed in 3] which circumvents the use of slow Monte-Carlo algorithm with 
relaxation time ^ L^. It relies upon the fact that the variables {i?i}o<j<L, given that (i?i+i — Hi) are weakly 
correlated random variables with their total sum equal to zero due to pbc, can be viewed as a Brownian bridge. To 
generate such a Brownian bridge, we first generate an ordinary random walk sequence Wi such that Wq = and 
Wj = Wi-i + rji, where rji^s are i.i.d. random variables drawn from Prob(?7 = x) ^ 9p{x) given in Eq. Next we 
express the Hi variables in terms of the Wi variables by the relation Hi = Wi — {i/L)yVL, such that they satisfy the 
same recursion relation as the Wi's, but this form ensures that the sequence of Hi's is periodic with period L and 
hence automatically enforces the pbc Hq = H^. One can show that this procdure yields the statistical weight for 
a configuration of Hi's, at least for large L (except for the p = 2 case where it is exact for all L). Once the actual 
height sequence Hi is generated, the relative height profile hi is then obtained from the relation hi = Hi — J^i Hi/L. 
For every sample of hi sequence, we measure the maximum hm and then obtain its histogram P(hm, L) from a total 
of 2 • 10^ configurations of the height profile. In addition, we also compute the width wl from w\ = J^ii^l) / 

In Fig. n we show a plot of ^/l2wj^ P{h„i, L) for p = 1 as a function of the rescaled variable hm/ \/^2w\ for 
different system sizes L — 256, 512 and 1024. The good collapse shows a good agreement with our analytical approach 
for p = \. Interestingly, our numerical simulations show that the leading behavior of P[h„i,L) as a function of 
hm/ \/l2wj^ (|41|l is in fact independent of p and given by the Airy-distribution function f{x) H39|l . This is depicted 
in Fig. [3 where we show a plot of yJVlwj^P{hm,L) for different values oi p — 0.5, 1, 1.5,2 and 3 and fixed L (for p 
fixed and different values of L one obtains similar results as for p = 1 in Fig. 

We have also investigated numerically the corrections to scaling. To do so, we introduce the quantity PncxtiKn, L) 

P„cxt(/J™, L) = ^/UwiPihm, L) - f{hrn/^/m^L) (43) 

For p = 1, according to our analytical calculations, the leading behavior of Pncxt(^mj^) is given by 
Pnc^tihm,L) - (2L)-i/2/'(/i„/(12w2)V2), 1^ pig. El we show a plot of {2LY/'^Pnc^t{hm,L) as a function of 
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hj(12w\)^ 



FIG. 1: \/l2'w\P{hm,L) as a function of fem/vl2w|^ 
for different system sizes L — 256, 512 and 1024 for 
p — 1. The solid line is the Airy-distribution function 
11391 1 evaluated using Mathematica. There is no fitting 
parameter. 




h„/(12w\)^ 



FIG. 2: ^Jl2wlP{hra,L) as a function of hm/^/Vh^ 
for p = 0.5, ...,3 and L — 512. The solid line is the 
Airy-distribution function 13911 evaluated using Math- 
ematica. There is no fitting parameter. 



the rescaled variable ft.m/\/l2w|_, for p = I, and different system sizes L = 64, 128 and 256. Here also, one obtains a 
good agreement with our analytical calculation H40|l . Finally, we have also computed numerically the corrections to 
scaling (|^ for different values of p. Our data (see Fig. ^ suggests that the leading correction to scaling, for any p, 
takes the following form 

Pnc.tihm,L) ^ ^i=/'(/i,„/./l2^) (44) 
where — 2. In Fig. 0] we show a plot of (/ipL)^/^Pnext(^m, as a function of the rescaled variable /im/\/l2i(;| 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 



h^/(12w\)"2 

FIG. 3: (2L)^''^Pncxt(/im,i) as a function of 
hm/ \/12w]^ for p — 1 and different system sizes L = 
64, 128, 256. The solid line is the derivative of the Airy- 
distribution function 1401 evaluated using Mathemat- 
ica. There is no fitting parameter. 
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h^/(12w\)^'2 

FIG. 4: {2Lf'''P^^^t{h^,L) as a function of 
hra.j \J 12iu| for different system sizes L — 128, 256 and 
for p — 1, 1.5, 2, with /ii = 2(exact), /xi.s = 3.8(1) and 
IJ.2 = 6.9(1). The solid line is the derivative of the 
Airy-distribution function 14011 evaluated using Math- 
ematica. 



for L — 128,256 and p — 1,1.5 and 2. In this figure, /ii.5 = 3.8(1) and /i2 = 6.9(1) are estimated to obtain the 
best collapse with f'{x). The good data collapse which we obtain is indeed in good agreement with the behavior 
in Eq. H44() . And therefore the corrections to scaling in Eq. (|41|l are also independent of p, up to a non universal 
amplitude ^p. 
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VI. CONCLUSION 

In this paper, we have studied the distribution of the maximal relative height /i„i in a class of one dimensional 
solid-on-solid lattice models defined by the Hamiltonian in Eq. (Q. We have provided a simple central limit argument 
to show that for arbitrary symmetric g{x) with a finite second moment cr^, the asymptotic MRH distribution P(hm, L) 
has a scaling form, P{hm,L) « (l/-\/T2wi)/ (/i^/vT^Wi) where the width w| w a'^L/12 for large L depends only 
on the variance ct^, but the scaling function /(x) is universal (independent of the details of the function g{x)) and is 
given by the Airy distibution function in Eq. H39f) . This argument is supported by exact calculation in the case of 
the Hamiltonian in Eq. ^ with p = 1 and also by numerical simulations of the model in Eq. (|SJ) for other values of 
p. Moreover, we have shown that even the leading correction to scaling function is also universal (up to an overall 
nonuniversal amplitude) and is given by the derivative of the Airy distribution function. 

There remain several open questions for future studies. In this paper we have studied symmetric g{x) functions 
with a finite second moment . We expect that the universality class of the asymptotic MRH distribution will change 
if is not finite. For example, if g{x) ~ a;"*^^"'"''-' for large x with < S < 2 that corresponds to the Levy walk of 
the interface, it would be interesting to compute the MRH distribution which we expect will depend on S. Another 
interesting case would be to study the MRH distribution with asymmetric g{x) that corresponds to interfaces with a 
drift. 

In this paper we have focused on very simple SOS Hamiltonians that have only 'kinetic' or 'interaction' part, but 
no onsite potential. There are important applications of the SOS Hamiltonians with an onsite attractive potential in 
wetting phenomena [2^ as well as in understanding the thermal denaturation of DNA molecules ■ In this case 
one considers a Hamiltonian of the form 

n = Y,[Fi\H^-H,+l\) + U{H,)] , (45) 

i 

where the heights Hi £] — oo,-t-oo[, F{x) — ~\n[g{x)] as in Eq. {TJ and U{Hi) represents a finite attractive onsite 
potential. One can ask: what is the MRH distribution in these models? In this case, both the heights iJ^'s as well 
as the relative heights hiS will be bounded variables in equilibrium. This follows from the transfer matrix analysis 
in section IV and the fact that any localized potential (far from the boundary) in one dimension will always give 
rise to a bound state. Consequently, the correlation function between the heights hiS will decay exponentially with 
distance. Thus, the MRH distribution in this case corresponds to calculating the distribution of the maximum of a 
set of 'weakly' correlated random variables. By considering 'blocks' of size of the correlation length, it is then easy 
to see that asymptotically (when the system size L is much larger than the corelation length) the MRH distribution 
will follow [l3| , after a suitable rescaling, the Gumbel law of the extreme of uncorrelated random variables |13| . 

On the other hand, if one considers the heights HiS to be positive variables. Hi > 0, i.e. one is considering a 
semi-infinite system, then the presence of a localized potential near the boundary H — leads to a roughening phase 
transition at a critical temperature Tc "2^. Below Tc, the interface is bound, i.e. localized near the surface at 
H = (with a finite L independent width for large L), but above Tc the system becomes unbound, i.e. rough with 
width wl ~ L^^^. Consequently, we would expect that the asymptotic MRH distribution will have a Gumbel form 
below Tc, but for T > Tc it will change to the Airy distribution function described in this paper. In short, the bound 
phase will be characterized by the Gumbel form and the rough phase by the Airy distribution function. It would 
be interesting to compute how this transition in the MRH distribution takes place across the critical temperature. 
Finally it would be interesting to compute the MRH distribution in the lattice SOS models in higher dimensions. 
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APPENDIX A: EXACT CALCULATION OF THE WIDTH wl IN THE SOS MODEL 

We wish to compute the width wl defined via 

^-lilihl) (Al) 

i=l 

where /i^'s are the relative heights 

h,^H,^^Y.^- (A2) 

i 

and HiS are the absolute heights in the SOS model defined by Eq. {Tj). Substituting Eq. (|A2|) in (jAl|l one gets 

- ^ E<^') - h H^H.H^)- (A3) 

i—l i<j 

The absolute height profile Hi's satisfy the recursion relation 

if, =i/,_i+6. (A4) 

It follows from Eq. jT)) that the set of random numbers {^i} with i = 1, 2, . . . , L are distributed according to the joint 
law 

Pm]^NLg{^i)cM2) ...gi^L)s(j2i}j (A5) 
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where N]^ is such that the joint distribution is normalized and g{x) is a symmetric, normahzed pdf with a finite 
second moment cr^, but otherwise arbitrary. The delta function ensures that the path is periodic, iJp = H]^. Thus, 
the variables ^^'s are only 'weakly' correlated due to the global constraint = 0. Since we are only interested in 

relative heights in Eq. (|A2p . we can, without any loss of generality, set iJo = Hl = 0. Thus, from Eq. IIA4p we have 

i 

H^^Y^^k. (A6) 

To calculate the correlation function {HiHj) as required in Eq. (|A3|) we need to compute the correlation function 
i^j^k) of the noise variables. This is easy to do. From the constraint = 0, it follows that (f^ ^i) — for all 

k. Using the isotropic property {i.e. (C/cCj) is the same for any pair of (fc, j) as long as fc ^ j), it follows that 



(aO) = -f^ (A7) 



where the onsite variance (^?) is independent of i and will be computed later. Eqs. I|A7|I and l|A6|l give the required 
height correlations 



{H,H,) = -^i(L-j) for i<j 



L- 1 

{&) 
L- 1 



jiL-i) for t>j. (A8) 



We now substitute this result in Eq. (|A3(I . We need to use one identity, namely, J2i<j ~ j) — + 1)(^ ^ — 
l)/24. Subsequently, a straightforward algebra gives 

= ^ i^f) ■ (A9) 

The only remaining task is to compute the onsite variance (^,^) where the set {Ci}'s is drawn from the joint 
distribution in Eq. IjASp . Let Pl{x) — Prob(^i = x) be the 'single site' probability density of the variable ^i, obtained 
by keeping £,i — x fixed and integrating the joint distribution in Eq. (jA5() over the rest {L — 1) variables. Evidently, 
it follows from the isotropic property that Pl{x) is independent of i. By integrating over the (L — 1) variables in Eq. 
HA5|) one gets 

Pl{x) = NLg{x)VL-ii-x) (AlO) 

where Vn{x) is the mutiple integral 



Vn{x) = J g{xi)g{x2) .■.g{xn)S 

g{xi)Vn-i{x - xi)dxi . (All) 



The recursion relation satisfied by Vn{x) above starts with the initial condition Vo{x) ~ 6{x) and has a simple 
convolution form. Hence it can be solved exactly by the Fourier transform method. Let, Vn{k) = Vn{x)e'^^^dx 
and g{k) — g{x)e'^''^dx. Taking the Fourier transform of Eq. (|Alip gives 

Mk) = g{k) Vn-i{k) = [g{k)r (A12) 
where we have used Vo{k) — 1. We invert the Fourier transform and substitute in Eq. I|A10|1 to get 

PL{x) = NLg{x) / [g{k)f~'e~'''^dk. (A13) 



The constant Nl is fixed by demanding PL{x)dx = 1. Using the symmety g{k) — g{^k) wc then get 

f°° \a{k)]'"^^ e^^^^ dk 
Pl{x) = g{x) ■'-"^ I: T . (A14) 

jZimfdk 
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We are now ready to compute the onsite variance, (^f) = J PL{x)x'^dx. We get from Eq. IIA14p 

where g"{k) = cPg/dk'^. The resuh m Eq. l|Il5|) can be further simphed via integration by parts. Substituting this 
result in Eq. (|A9|I we get our final exact result, valid for all L, 

where g'{k) — dg/dk. 

Let us now consider a few special cases. For the Gaussian case, i.e. for p = 2 in Eq. ®, we have g{k) — 
exp[-cr2^2/2]. From Eq. (|IT6)l we get for all L > 1 

For the p = 1 case in Eq. (0, we have g{x) = K exp[~K\x\]/2. The Fourier transform gives g{k) = K"^ /{K"^ + fc^). 
Thus in this case, = —g"{0) = 2/K'^. Substituting in Eq. (|A16|) and performing the integral we get, for all L > 1 

wr = P -^o- . (A18) 

^ 24L(i + l) ^ ' 

In general for arbitrary symmetric g{x) with a finite second moment cr^, one gets from Eq. (|A16|I the following 
large L behavior, 

wl^a'^. (A19) 

This can be seen from the fact that for large L, the Fourier integrals in Eq. (jA16p will be dominated by the 
contributions from the small k limit. For small k, g{k) « 1 — a'^k^/2 + ... « exp[— (T^fc^/2]. Substituting this result in 
Eq. (jA16|l and performing the integrals one arrives at the asymptotic behavior in Eq. I|A19|) . 



